Under consideration for publication in J. Fluid Mech. 



1 



A theory for the emergence of coherent 
structures in beta-plane turbulence 

Nikolaos A. Bakas^f, and P. J. loannou^ 

^National and Kapodistrian University of Athens, Build. IV, office 34, Panepistimiopolis, 

Zografos, Athens, Greece 

(Received 22 March 20f3; revised ?; accepted ?.) 

Turbulent flows are observed to self-organize into large scale structures such as zonal 
jets and coherent vortices. The simplest model that retains the relevant dynamics is a 
barotropic Row in a /3-plane channel with turbulence sustained by random stirring. Non- 
linear integrations of this model show that as the energy input rate of the forcing is 
increased, the homogeneity of the flow is broken first by the emergence of non-zonal, co- 
herent, westward propagating structures and second for larger energy input rates by the 
emergence of zonal jets. We use a non-equilibrium statistical theory, Stochastic Struc- 
tural Stability Theory (SSST), to address the emergence of these coherent structures. 
Employing the tools of SSST, we construct a coupled dynamical system governing the 
joint evolution of the coherent flow and of the second order statistics of the turbulent 
eddies. We then treat the structural stability of a homogeneous equilibrium with no mean 
flow analytically. We find that non-zonal and zonal coherent structures appear as the re- 
sult of structural instability and equilibrate at finite amplitude. We also find that SSST 
accurately predicts the critical energy input rates for the emergence of both non-zonal 
and zonal coherent structures in the non-linear integrations as well as the characteristics 
(scale, amplitude and phase speed) of these structures. 

1. Introduction 

Atmospheric and oceanic turbulence is commonly observed to be organized into spa- 
tially and temporally coherent structures such as zonal jets and vortex rings. Examples 
from planetary turbulence include the banded jets and the Great Red Spot in the Jo- 
vian atmosphere (Ingersoll 1990; Vasavada & Showman 2005), as well as the latent jets 
in the Earth's ocean basins (Maximenko et al. 2005) and the ocean rings shed by the 
meandering of the Gulf-Stream in the western Atlantic Ocean (Chelton et al. 2007). Lab- 
oratory experiments and numerical simulations of both decaying and forced turbulence 
have shown that these coherent structures can spontaneously appear out of a background 
of homogeneous turbulence and persist for a very long time despite the presence of eddy 
mixing (VaUis & Maltrud 1993; Cho & Polvani 1996; Weeks et al. 1997; Read et al. 2004; 
Espa et al. 2010). 

The simplest model that captures the turbulent dynamics and its interaction with the 
coherent structures, is the stochastically forced barotropic vorticity equation on the sur- 
face of a rotating planet or on a /3-plane (a plane tangent to the surface of a rotating 
planet in which differential rotation is taken into account). A large number of numerical 
simulations of this model have shown that robust, large scale zonal jets emerge in the 
flow and are sustained at finite amplitude (Williams 1978; Vallis & Maltrud 1993; Nadiga 
2006; Danilov & Gurarie 2004; Galperin et al. 2006). In addition, large scale westward 
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propagating coherent waves were found to coexist with the zonal jets (Sukariansky et al. 
2008; Galperin et al. 2010). These waves were found to cither obey a Rossby wave disper- 
sion, or propagate with significantly lower phase speeds. The slowly propagating waves 
typically have low zonal wavenumbers and were found in a parameter regime in which 
strong, robust jets dominate. These waves that were termed as satellite modes (Danilov 
& Gurarie 2004) or zonons (Sukariansky et al. 2008), appear to be sustained by non- 
linear interactions between Rossby waves. However the mechanism for their excitation 
and maintenance remains elusive. The goal in this work is to develop a non-equilibrium 
statistical theory that can predict the emergence of both zonal jets and non-zonal coher- 
ent structures and can capture their characteristics. 

The tendency for formation of large scale structures in planetary turbulence is a conse- 
quence of the approximate energy and vorticity conservation in two dimensional or quasi 
two dimensional flows that implies an energy transfer from small to large scales (Fjortoft 
1953). Rhines (1975) found using the barotropic model that the non-linear eddy-eddy in- 
teractions that arc local in wavenumber space, lead to an isotropic inverse energy cascade 
similar to two-dimensional Navier-Stokcs turbulence (Kraichnan 1967). However, Rhines 
(1975) found that the isotropic cascade is inhibited in the region in wavenumber space 
in which weakly interacting Rossby- waves dominate. As a result it is anisotropized and 
continues through a narrow region in wavenumber space, transferring energy to zonal 
jets (Vallis & Maltrud 1993; Nazarenko & Quinn 2009). The energy cascade is flnally 
arrested at a meridional scale that is dictated by friction (Smith et al. 2002; Sukariansky 
et al. 2007). However, observations of the atmospheric midlatitude jet (Shepherd 1987) 
and numerical analysis of simulations (Nozawa & Yoden 1997; Huang & Robinson 1998; 
Huang et al. 2001), showed that the large scale jets are maintained through spectrally 
non-local interactions rather than by a local in wavenumber space cascade. Moreover, the 
scales and characteristics of the non-zonal coherent structures that also emerge cannot 
be explained by the phenomenological description of the inverse turbulent cascade. As 
a result, a comprehensive theory for the emergence of coherent structures needs to be 
developed. 

Since organization of turbulence into coherent structures involves complex non-linear 
interactions among a large number of degrees of freedom, an alternative approach for 
gaining an understanding for the tendency towards self-organization of turbulent flows 
is to use statistical mechanics, an approach pioneered by Miller (1990) and Robert & 
Sommeria (1991) in what is now known as Robert-Sommcria-Miller (RSM) theory. The 
RSM theory builds upon the work of Onsager (1949) that explains self-organization of 
turbulence in terms of the equilibrium statistical mechanics of a set of point vortices. 
The main idea is to use equilibrium statistical mechanics to flnd a solution of the ideal 
and unforced Euler equations that maximizes a proper measure of entropy under the 
restrictions imposed by all conserved quantities in the equations. The coherent structures 
that emerge out of this statistical analysis for two dimensional and quasi-geostrophic flows 
are either large scale vortices (Chavanis & Sommeria 1998) or jets (Bouchet & Sommeria 
2002; Venaille & Bouchet 2011) (see also a recent review by Bouchet & Venaille (2012)). 
A recent study has also shown correspondence of the theoretical predictions with non- 
linear simulations in the limit of weak forcing and dissipation (Bouchet & Simonnet 
2009). However, the relevance of these results in planetary flows that are strongly forced 
and dissipated and are therefore out of equilibrium remains to be shown. 

A non-equilibrium statistical theory, that is termed as Stochastic Structural Stability 
Theory (SSST) (Farrell & loannou 2003) or Second Order Cumulant Expansion theory 
(CE2) (Marston et al. 2008) and can treat turbulent flows out of equilibrium was de- 
veloped during the last decade and was applied to macroscale barotropic and baroclinic 
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turbulence in planetary atmospheres, plasmas and astrophysical flows (Farrell & loan- 
nou 2003, 2007, 2008; Marston et al. 2008; Farrell & loannou 2009a, 6, c; Marston 2010; 
Tobias et al. 2011; Srinivasan & Young 2012; Marston 2012). This theory is based on 
two building blocks. The first is to decompose the dynamical variables into the sum of a 
mean value that represents the coherent flow and fluctuations that represent the turbu- 
lent eddies. The cumulants of the dynamical variables containing the information on the 
mean values of the variables and their correlations are then formed and the equations 
governing the evolution of the cumulants are obtained. These equations govern the joint 
evolution of the mean coherent flow (first cumulant) and the ensemble eddy statistics 
(higher order cumulants) . The second building block is to either parameterize the terms 
involving the third cumulant as stochastic excitation and enhanced dissipation (Farrell 
& loannou 1993a, 6, c; DelSole & Farrell 1996; DelSole 2004) or set the third cumulant 
to zero (Marston et al. 2008; Tobias et al. 2011; Srinivasan & Young 2012). Restriction 
of the dynamics to the first two cumulants is equivalent to neglecting the eddy-eddy 
interactions in the fully non-linear dynamics and retaining only the interaction between 
the eddies with the instantaneous mean flow. This closure results in a non-linear au- 
tonomous dynamical system that governs the evolution of the coherent mean flow and 
the second order eddy statistics. A similar approach was also followed by DubruUe and 
collaborators (DubruUe & Nazarenko 1997; Laval et al. 2003) to describe the interac- 
tion of coherent vortical structures with turbulence. This second order closure has been 
shown to produce accurate quadratic eddy statistics and mean flows in studies of plan- 
etary turbulence (DelSole & Farrell 1996; DelSole 2004; O'Gorman & Schneider 2007). 
The reason for the accurate predictions of this closure is that it captures the non-local 
in wavenumber space interaction between the perturbations and the large scale structure 
that is found to sustain the mean flows in planetary turbulence. 

The averaging operator needed to separate the mean from the fluctuations in this the- 
ory can be defined differently, depending on the physical situation addressed. In most 
earlier studies of SSST which addressed the interaction of zonal jets with turbulence, 
a zonal average was employed and the emergence and equilibration of turbulent jets 
in barotropic and baroclinic rotating atmospheres was investigated (Farrell & loannou 
2007, 2008, 2009a). Direct comparisons of the predictions of SSST using the zonal mean 
interpretation for the emergence and equilibration of jets with fully non-linear barotropic 
models were recently reported (Srinivasan & Young 2012; Constantinou et al. 2013; To- 
bias & Marston 2013). It was found that SSST can accurately predict the structure of 
zonal flows in the non-linear simulations. However, with this interpretation, the non- 
zonal waves are treated as incoherent and their emergence and characteristics cannot be 
studied. 

In order to investigate the dynamics of the coherent non-zonal structures, we adopt in 
this work the more general interpretation that the ensemble average represents a Reynolds 
average with the ensemble mean representing coarse-graining, an interpretation that has 
also been recently adopted in SSST studies of baroclinic turbulence (Bernstein 2009; 
Bernstein & Farrell 2010). With this interpretation of the ensemble mean, we obtain the 
statistical dynamics of the interaction of both zonal and non-zonal coherent structures 
with stochastically forced turbulence on a barotropic /3-planc channel, with the goal of 
addressing their emergence and characteristics. We find that the turbulent equilibrium 
that is homogeneous, is structurally unstable when the energy input rate is above a 
threshold and both zonal and non-zonal coherent structures emerge. We also show that 
the characteristics of these structures observed in the non- linear simulations are predicted 
by SSST. 

This paper is organized as follows. In section 2 we present the characteristics of the 
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zonal and non-zonal coherent structures that emerge in non-linear simulations of the 
turbulent flow. In section 3 we derive the SSST system that governs the evolution of the 
ensemble mean coherent structures (first cumulant) and the associated eddy statistics 
(second cumulant). In section 4 we study the instability of the corresponding homoge- 
neous equilibrium, analyzing the unstable structures and their dispersion relation and 
we investigate the equilibration of the instabilities in section 5. The predictions of SSST 
are then compared to the results of the non-linear simulations in section 6 and we finally 
end with a brief discussion of the obtained results and our conclusions in section 7. 



2. The emergence of coherent structures in non-Unear simulations of 
a barotropic flow 

Consider a nondivergent barotropic flow on a /3-plane with cartesian coordinates x = 
(x,y). The velocity field, u = {u,v), is given by {u,v) = {—ipy,ipx), where ip is the 
streamfunction. Relative vorticity (^{x,y,t) ~ Aip, evolves according to the non-linear 
equation: 

{dt+u-V)C + Pv = ~rC-i^A''C + r, (2.1) 
where A = d^^ + dyy is the horizontal Laplacian, /3 is the gradient of planetary vorticity, 
r is the coefficient of linear dissipation that typically parameterizes Ekman drag and v is 
the coefficient of hyper-diffusion that dissipates the energy flowing into unresolved scales. 
The forcing term is necessary to sustain turbulence and serves as a parameterization of 
processes that are missing from the barotropic dynamics, such as small scale convection 
or baroclinic instability. We will consider the flow to be on a doubly periodic channel of 
size 27r x 27r. 

As in many previous studies the exogenous excitation will be assumed to be a 
temporally delta correlated and spatially homogeneous and isotropic random stirring 
with a two-point, two-time correlation function of the form: 

{.r{xi,yi,ti)f''{x2,y2,t2)) = S{t2 - ti)E{xi,X2,yi,y2)- (2.2) 

The brackets denote an ensemble average over the different realizations of the forcing 
and the spatially homogeneous covariance, S, in the doubly periodic channel is: 



S(xi,X2,2/i,y2) - ^^S(fc,Z)e"=^— (2.3) 

k I 

with the X, y wavenumbers, k and I, taking all integer values. The Fourier amplitude 



eKf f 1, for Wk^ + P-Kf \ i^AKf 



" AKf [ 0, for \V^TP -Kf\>AKf ' ^^'^^ 

is chosen so that the excitation injects energy at rate e in a narrow ring in wavenumber 
space with radius Kf and width AKf. 

Equation (2.1) is solved using a pseudospectral code with a 128 x 128 resolution and 
a fourth order Runge-Kutta scheme for time stepping. The parameters used are (3 = 10, 
r ~ 0.01, u = 1.19 • 10~^, Kf = 10 and AKf = 1, and are typical of planetary flows. 
In addition, the results presented were found to be insensitive on the chosen parameters 
and the forcing protocol. The non-linear system reaches a statistical equilibrium at about 
t — 10/r. Following previous studies (Galperin et al. 2006), the integration was carried 
until t = lOO/r in order to collect accurate statistics and the last 80/r time units were 
used for calculating the time averages. To illustrate some of the characteristics of the 
turbulent flow and the emergence of structure, we consider two indices that measure the 
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power which is concentrated at scales larger than the scales forced. The first is the zonal 
mean flow index defined as in Srinivasan & Young (2012), as the ratio of the energy of 
zonal jets over the total energy 

where E(k, I) is the time averaged energy power spectrum of the flow at wavenumbcrs 
(fc, I). The second is the non-zonal mean flow index defined as the ratio of the energy of 
the non-zonal modes with scales lower than the scale of the forcing over the total energy: 

y:ki:K<K,-AK, m o - Ei Hk = o, /) 

nzmf = —. (2.6) 

Em mi) 

If the structures that emerge are coherent, then these indices quantify their amplitude. 
Figure 1 shows both indices as a function of the energy input rate e. Remarkably, both 
indices exhibit sharp increases at critical energy input rates, indicating the occurrence 
of regime transitions in the flow. For e smaller than a critical value ec, the turbulent 
flow is homogeneous and remains translationally invariant in both directions and both 
indices are zero. When e > ec, non-zonal structures that have scales larger than the scale 
of the forcing form, as indicated by the increase in the nzmf index. The time averaged 
power spectrum shown in figure 2(a) for e = 4ec, is anisotropic with a pronounced peak 
at (|fc|, |/|) = (1,5). This peak corresponds to a structure with the corresponding scale 
that is evident in the vorticity field evolution. This is illustrated by the appearance of 
a structure with (|fc|, |/|) = (1,5) in the snapshot of the streamfunction field shown in 
figure 3(a). The Hovmoller diagram of flj{x,y = 7r/4,i) plotted in figure 3(b) shows that 
this structure is coherent and propagates in the retrograde direction. The slope in the 
diagram corresponds to the phase speed of the waves, which is found to be approximately 
the Rossby wave phase speed for (fc, I) = (1,5). This structure remains dominant for larger 
input rates as well, with its energy comprising over 60 % of the total power in non-zonal 
structures for e/ec < 15. Therefore the increase in the nzmf index observed in figure 1 
signifies the emergence of coherent structures that break the translational symmetry of 
the turbulent state. 

For e > e„i the increase in the zmf index, shown in figure 1, indicates a second regime 
transition in the flow with the emergence of robust and coherent zonal jets. For e = 3.3e„i 
(i.e. e/cc = 50) the spectrum, shown in figure 2(b), has significant power at the zonal 
structures with (fc, |Z|) = (0,4). These peaks correspond to coherent zonal jets as il- 
lustrated by the Hovmoller diagram of the zonally averaged streamfunction shown in 
figure 3(c). However, there is significant power in non-zonal structures (in this case with 
wavenumbcrs (|A:|, |^|) = (1,4) and (|A:|, |/|) = (1,5)), a characteristic that is also revealed 
by the high values of the nzmf index for large energy input rates. The Hovmoller diagram 
in figure 3(d) reveals that the non-zonal structures are coherent and propagating in the 
retrograde direction. However, the phase speed calculated from the diagram is lower than 
the corresponding Rossby wave speed for both (|fc|, \l\) = (1, 4) and (|fc|, |/|) — (1, 5). Sim- 
ilar Rossby-like, westward propagating coherent structures were also reported recently 
in numerical simulations of the barotropic vorticity equation on the sphere (Sukariansky 
et al. 2008; Galperin et al. 2010). In agreement with the results presented in this work 
these large scale waves contain a significant amount of energy. In the regime in which 
zonal jets are absent or weak (termed as friction dominated in these studies) these waves 
were found to follow the Rossby wave dispersion. In the regime in which strong zonal 
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Figure 1. The zmf (red lines) and nzmf (blue lines) indices defined in (2.5) and (2.6) respectively, 
as a function of energy input rate e/ec for the non- linear (solid lines) and quasi- linear (dashed 
lines) integrations. The critical value ec is the energy input rate at which the SSST predicts 
structural instability of the homogeneous turbulent state. Zonal jets emerge for e > e„;, with 
Eni = 15ec for the parameters in the simulation. 



jets dominate the flow (termed as zonostrophic regime in these studies), the waves prop- 
agate with markedly lower phase speeds. These waves were therefore termed as linear 
Rossby waves in the former and as satellite modes (Danilov & Gurarie 2004) or zonons 
(Sukariansky et al. 2008) in the latter regime. 

We will show next that the emergence and characteristics of both the zonal and the 
non-zonal coherent structures can be accurately predicted by considering the stability of 
a particular second order closure of the turbulent dynamics. This second order closure 
results in a non-equilibrium statistical theory, called Stochastic Structural Stability The- 
ory (SSST) or Second Order Cumulant Expansion theory (CE2) (Farrell & loannou 2003, 
2007; Marston et al. 2008; Bakas & loannou 2011; Srinivasan & Young 2012; Marston 
2012), that addresses the emergence of structure in planetary turbulence. 
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Figure 2. Time averaged energy power spectra, log(_B(fe, I)), obtained from the non-linear sim- 
ulation of (2.1) at (a) e/tc = 4 and (b) e/ec = 50 or e = 3.3e„;. In (a) the flow is dominated by a 
(|fc|, |/|) = (1,5) non-zonal coherent structure. In (b) the flow is dominated by a coherent zonal 
flow at (fc, |;|) = (0,4). 



3. Formulation of Stochastic Structural Stability Theory 

SSST describes the statistical dynamics of the first two equal time cumulants of (2.1). 
The first cumulant is the ensemble mean of the vorticity Z(x, <) = (C). The second 
cumulant C(xi,X2,t) = {('1(2)^ is the two point correlation function of the vorticity 
deviation from the mean C^' = Q — Zi, where we have used the shorthand (^i = Q(xi,t), 
with i = 1, 2 to refer to the value of the relative vorticity at the points Xi = (xi,yi). 
In most earlier studies of SSST, the ensemble average was assumed to represent a zonal 
average. With this interpretation of the ensemble average, the non-zonal structures are 
treated as incoherent motions and the theory can only address the emergence of zonal 
jets. In order to address the emergence of coherent non-zonal structures in turbulence, 
we adopt in this work the more general interpretation that the ensemble average is a 
Reynolds average over the fast turbulent motions. This interpretation of the SSST has 
been adopted recently in studies of atmospheric blocking patterns in baroclinic two-layer 
turbulence by Bernstein (2009) and Bernstein & Farrell (2010). With this definition of the 
ensemble mean, we seek to obtain the statistical dynamics of the interaction of the coarse- 
grained ensemble average field, which can be zonal or non-zonal coherent structures 
represented by their vorticity Z, with the fine-grained incoherent field represented by its 
vorticity covariance C. 

The equations governing the evolution of the first two cumulants are obtained as fol- 
lows. Under the decomposition of vorticity into an ensemble mean and a deviation from 
the mean, (2.1) is split into two equations governing the evolution of the eddy (deviation 
from the mean) vorticity and the vorticity of the coherent structures Z: 

[dt + U- V) C + {13 + Zy)v' + u'Z, = -rC - u^'C + r + r\ (3.1) 

/ 



(a* + LT • V) Z + = -V • {u'O - rZ - i^A^Z, (3.2) 
where u' = [u', v'] — [—dy^p', dxip'] denotes the eddy velocity vector, U — [U, V] = 
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Figure 3. (a) Snapshot of the streamfunction tp{x,y,t) and (b) HovmoUer diagram of 
ip{x,y — 7r/4, obtained from non-linear simulation at e/tc = 4. The thick dashed line in 
(b) correspond to the phase speed obtained from the stability equation (4.2). (c) HovmoUer 
diagram of the a;-averaged tp(y,t), which is obtained from non-linear simulations at e/ec = 50. 
(d) HovmoUer diagram of ip{x,y = 7r/4, t) obtained from the non-linear simulation at e/ec = 50. 
The thick dashed line corresponds to the phase speed obtained from the stability equation (4.2). 
At y = 7r/4 at which tp is calculated, the jet is maximum. The non-zonal structures remain 
coherent at other locations as well. 



[-~dy'9 , dx"^] is the corresponding ensemble mean velocity vector, 

= (it' • vc') - It' • vc', 



(3.3) 



is the forcing term from the non-linear interactions among the turbulent eddies and 
f = + represents the total eddy forcing. The ensemble average vorticity fluxes 
(u'C) can be expressed in terms of the second cumulant of vorticity as: 



(3.4) 



where the variables ai and the operators in (3.4) are evaluated at points x^, the operator 
inverts vorticity into the streamfunction field {ipi = ^T^C'i) ^^'^ the subscript Xi = 
X2 means that the expression in parenthesis is calculated at the same point. As a result, 
the first cumulant evolves as: 

dtZ + UZx + V{l3 + Zy) = dx {dy,A^'C)^^^^^-dy {O, , ' c) ^^^^^ - T Z - Z . (3.5) 

The equation governing the evolution of the second cumulant C can be obtained by first 
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rewriting (3.1) in the compact form: 

dtC - AC + f^: (3.6) 

where 

A^ = -UA^ - V,dy^ -{13 + Zy,)d^^Ar^ + Z^^dyAl^ (3-7) 
governs the dynamics of hnear perturbations about the instantaneous mean flow XJ . 
Multiplying (3.6) for 9tCi by C2 ^"^^ (3-6) for dtC^^ by Ci; adding the two equations and 
taking the ensemble average yields: 

dtC = (^1 + A^)C + (AC^ + , (3.8) 

where 

(/ic^ + /2C0 = (/rc^ + /2C0 + (/r'c^ + /2"'co 
= (/rc2+/2 + 

+ - ^^..3) + [(C.3 - 9lyJ ^2'^]^..=., ' (3-9) 

and r = (C1C2C3) is the third cumulant. The first term on the right hand side of (3.9) 
is the correlation of the external forcing with vorticity, while the other two terms 
involve the third cumulant that describes the eddy-eddy interactions. Previous studies 
addressing the interaction of turbulent eddies with zonal jets in baroclinic turbulence, 
as well as the interaction of coherent vortices with small scale turbulence, have shown 
that several important features of the coherent flow as well as accurate eddy statistics 
are obtained by either neglecting or suitably parameterizing the eddy-eddy non-linearity 
as stochastic forcing and enhanced dissipation (Farrell & loannou 1993a; DelSole & 
FarreU 1996; DubruUe & Nazarenko 1997; Laval et al. 2000; DelSole 2004; O'Gorman & 
Schneider 2007; Marston et al. 2008). This is equivalent to setting the third cumulant 
to zero, or parameterizing the last two terms in (3.9) as stochastic forcing with a known 
correlation function. In this work we will neglect the third-order cumulants and show that 
the SSST theory with this approximation can accurately predict the emergence of large 
scale structures in the flow. We therefore assume that / is the delta correlated external 
forcing, f^, with (/1C2 + /2C1) — (/1/2) = With this approximation the second order 
statistics evolve according to: 

dtC = {Ai+A2)C + ^. (3.10) 

Equations (3.5) and (3.10) form a closed deterministic system that governs the joint 
evolution of the coherent flow field and of the second order turbulent eddy statistics. 
This second order closure is the basis of Stochastic Structural Stability Theory (Farrell 
& loannou 2003). The SSST system can have fixed points, limit cycles or chaotic attrac- 
tors. Examples of the attractor of this system can be found in the SSST description of 
the organization of geophysical and plasma turbulence into zonal jets (Farrell & loannou 
2003, 2008, 20096), as well as in the SSST description of blocking patterns in the at- 
mosphere (Bernstein & Farrell 2010). The fixed points and C^, if they exist, define 
statistical equilibria of the coherent structures with vorticity, Z^, in the presence of an 
eddy field with covariance, . The structural stability of these turbulent equilibria that 
can be investigated in SSST, addresses the parameters in the physical system which can 
lead to abrupt reorganization of the turbulent fiow. Specifically, when an equilibrium of 
the SSST equations becomes unstable as a physical parameter changes, the turbulent flow 
bifurcates to a different attractor. In this work, we show that coherent structures emerge 
as unstable modes of the SSST system and equilibrate at finite amplitude. The predic- 



10 



N. A. Bakas and P. J. loannou 



tions of the SSST system regarding the emergence and characteristics of the coherent 
structures are then compared to the non-hnear simulations. 



4. SSST instability and emergence of finite amplitude large scale 
structure 

It can be seen that the homogeneous equihbrium with no mean flow 

Z^ = 0, C"^^, (4.1) 

is a fixed point of the SSST system (3.5) and (3.10) in the absence of hyperdiffusion 
(see Appendix A) . The stability of this homogeneous equilibrium, can be addressed by 
performing eigenanalysis of the SSST system linearized about the equilibrium. Because of 
the absence of coherent mean flow and the homogeneity of we can seek eigensolutions 
in the modal form 5Z = Z„„,e'"^+™'!'e'"* and 5C C„™(x, y)e'"^+'™«e'"*, where x = 
Xi — X2, X = {xi + a;2)/2, y = yi — y2, y ~ [yi + y2)/2, n and m are the x and y 
wavenumbers of the eigenfunction and cr = (7^ + if j is the eigenvalue with ct^ = Re(cr), 
Gi = Im((T) being the growth rate and frequency of the mode respectively. The eigenvalue 
a satisfies the equation: 

(mfc - nl) [nm{kl - iD + (m^ - n^)k+l+] K^{K^ - N^)S^{k, I) 
, 2if3k+{k+n + l+m) - in/3 {K"^ + A'2) /2 + (cr + 2r)K^K^ 
= TT{a + r)N^ -mn/3, (4.2) 

where 

-1 /'OO /'OO 

^''(A;,/) = ^^^y J S(5,y)e-'=*""Midy , (4.3) 

is the Fourier transform of the streamfunction covariance S = (V'iV'2) equilibrium, 
K"^ = k"^ + Kl = {k + n)2 + {l + m)2, = n"^ + m^, k+ = k + n/2 and 1+ = 
I + m/2 (cf. Appendix A). For zonally homogeneous perturbations with n = 0, (4.2) 
reduces to the eigenvalue relation derived by Srinivasan & Young (2012) for the emergence 
of jets in a barotropic /3-plane. Eigenvalue relation (4.2) was derived for a flow that 
extends to infinity. For the periodic channel considered in the non-linear simulations, 
the corresponding eigenvalue relation is readily obtained by substituting the integrals in 
(4.2) and (4.3) with summation over integer values of k and I. We non-dimensionalize 
the eigenvalue relation using the dissipation time scale 1 /r and a typical forcing length 
scale Lf and rewrite (4.2) in the general form: 



' (n,m) 



.g(/3,e). (4.4) 



For a given spectral distribution of the forcing, (4.4) gives the eigenvalue a = a/r for 
each wavenumbers {h,rh) = Lf(n,m) as a function of the planetary vorticity gradient 
/3 = (3Lf/r and the energy injection rate e = e/{r^L'j). 

We consider the case of a ring forcing that injects energy at rate e at the total wavenum- 
ber Kf. 



S(fc, /) = 2eKfS{vW+P - Kf), (4.5) 

that is an idealization of the ring forcing (2.4) considered in the non-linear simulations 
and solve (4.4) numerically. For small values of the energy input rate, the growth rate 
CTr is negative for all structures and the homogeneous equilibrium is stable. At a crit- 
ical ic the homogeneous flow becomes SSST unstable, symmetry breaking occurs and 
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exponentially growing coherent structures emerge. The critical value, Ec, is calculated by 
first determining the energy input rate et{n, rh) that renders wavenumbers (n, m) neutral 
{^r{n,m) = O) j Stnd then by finding the minimum energy input rate over all wavenumbers: 
ic = min(fi ,^)et. This minimum input energy £c as a function of /3 is shown in figure 4. 
The absolute minimum energy input rate required is ic = 67 and occurs at /3mm = 3.5. 
For /3 ^ Pmim the structures that first become marginally stable are zonal jets (with 
n = 0). The critical input rate increases as Cc ^ for /3 — > (in agreement with the 
findings of Srinivasan & Young (2012)) and the homogeneous equilibrium is structurally 
stable for all excitation amplitudes when /3 = 0. The structural stability for /? = is an 
artifact of the assumed isotropy of the excitation and in the presence of anisotropy the 
critical input rate, ec, saturates to a finite value as /3 — )• (Bakas & loannou 2011, 2013). 
For j3 > Pmin, the marginally stable structures are non-zonal and ic grows as ^ 0^1"^ 
for /3 —s- OD. Since the critical forcing for the emergence of zonal jets (also shown in figure 
4), increases as tc ^ 0^ for /3 — >■ cjo (Srinivasan & Young 2012), for large values of /3 non- 
zonal structures first emerge and only at significantly higher e zonal jets are expected to 
appear. Investigation of these results with other forcing distributions revealed that these 
results are independent of the isotropy of the forcing. Contours of the maximum growth 
rate of the SSST instability, (imax = 'Cii^'y^(n,m)^r are also shown in figure 4 as a function 
of (e, /3). For a given /?, the maximum growth rate increases monotonically with larger 
energy input rates, while for a given level of excitation e„i the maximum growth rate 
occurs for a finite that satisfies roughly tm ~ 30/3^ (represented by the thick dotted 
line in the figure). 

For e > Ec the unstable mode eigenvalues satisfy the relations: 

^(-h,m) = 0'ln,m)^ ^^'^ <^{h,-m} = ^{ii,m), (4.6) 

implying that the growth rates depend on |?t,| and |m| (the symmetries in (4.6) are 
satisfied for all cases treated in this paper, but arc not in general valid). As a result, the 
plane wave 6Z = cos{nx + my) and an array of localized vortices SZ = cos(na;) cos(my), 
have the same growth rate, despite their different structure. 

We first consider the case /3 — 1, e = 2ec, corresponding to the point marked as 
5a in the (e, /3) regime diagram shown in figure 4. The growth rate of the maximally 
growing eigenvalue, ar, and its associated frequency of the mode, (Ti, arc plotted in 
figure 5a as a function of |n| and |77i|. We observe that the region in wavenumber space 
defined roughly by < |n| < 1/2, and 1/2 < |m| < 1 is unstable, with the maximum 
growth rate occurring for zonal structures (n = 0) with |m| ~ 0.8. The frequency of the 
unstable modes is in general non-negative (ct^ ^ 0) and is equal to zero only for zonal jet 
perturbations (n = 0). Using the symmetries (4.6), this implies that real unstable mean 
flow perturbations SZ propagate in the retrograde direction ii n ^ and are stationary 
when n = 0. As e increases the instability region expands roughly covering the sector 
1/2 < |A^| < 1 and a second instability branch with smaller growth rates appears for 
\N\ > 1. This is illustrated in figure 5b showing the growth rate for e = lOcc (marked 
as 5b in figure 4). Note also that for both values of the energy input rate, the zonal 
structures have a larger growth rate than the non-zonal structures, a result that holds 
for any e when /S < (3min- 

For /3 > Pmin the non-zonal structures have always larger growth rate. This is il- 
lustrated in figures 6 and 7, showing the growth rates and frequencies of the unstable 
modes for /3 = 10 and /3 = 100 respectively. For larger /3 values there is tendency for the 
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Figure 4. The critical energy input rate Ic for structural instability (thick solid line) and 
the critical energy input rate for structural instability of zonal jets (solid line) as a function 
of /3. The behavior of these critical values for large and small /? is indicated with the dashed 
asymptotes e = 23/3~^ for /3 ^ 1, e = 11^^''^ and e = 0.5/3^ for the emergence of non-zonal and 
zonal structures respectively for /3 ^ 1. Above the thick (thin) line non- zonal (zonal) coherent 
structures emerge. The thin dotted vertical line /3 = Pmin separates the unstable region: for 
/S < /3min the zonal structures grow the most, whereas for /? > j3min the non-zonal structures 
grow the most. Also shown are the contours (thick dashed lines) of the maximum growth rate 
d'max (with contour values corresponding to log{amax)) ■ The thick dotted line e — 30/3'^ is the 
locus of the points on which the maximum ct^ occurs for each e. The crosses indicate the e and 
P values for which the dispersion relation of the unstable modes is shown in figures 5-7. 



frequency of the unstable modes to conform to the corresponding Rossby wave frequency 

an = .rj^ , (4.7) 

a tendency that docs not occur for smaller A comparison between the frequency of the 
unstable mode and the Rossby wave frequency is shown in figures 7(c), (d) in a plot of 
di/ffR. For slightly supercritical e, the ratio is close to one and the unstable modes satisfy 
the Rossby wave dispersion relation. At higher supercriticalities though, ct^ departs from 
the Rossby wave frequency (by as much as 40% for the case of e = SOcc shown in figure 
7(d)). 
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(a) 




Figure 5. Dispersion relation of the unstable modes for 13 = 1. The energy input rate is (a) 
e = 2ec and (b) e = lOtc- The contours show the growth rate ar and the shading shows the 
frequency ct, of the unstable modes. 



5. Equilibration of the SSST instabilities 

Wc now demonstrate the equilibration of the instabilities by the discretizcd SSST sys- 
tem (3.5), (3.10) in a doubly periodic channel of size 27r x 2tt. We consider the parameter 
values chosen in the non-linear simulations discussed in section 2 (/3 = 10, r = 0.01, 
1/ = 1.19 • 10~^, Kf = 10 and A/v/ ~ 1). For these parameters, /? = 100 and therefore 
the integration is in the parameter region of figure 4 in which the non-zonal structures are 
more unstable than the zonal jets. We first consider the energy input rate e = 4ec which 
corresponds to the first case presented in section 2 and is marked as 8a in figure 4f . The 
growth rates of the coherent structures for integer values of the wavenumbers. n and m 
are calculated from the discrete version of equation (4.2), because of the 2tt periodicity 
of the channel, with the addition of a hyperdiffusive dissipation term in equation (4.2). 
The resulting growth rates for this energy input rate are shown in figure 8(a). For these 
parameters the zonal jet perturbations are stable and are not expected to emerge, while 
a large number of non-zonal modes are unstable with the perturbation (n,m) = (1,5) 
growing the most. At t = 0, wc introduce a small random perturbation, whose stream- 
function is shown in figure 9(a). After about t — 40/cr(i 5), where cr(i.5) is the growth rate 
of (ri, m) — (1, 5), a checkerboard perturbation of the form Z = cos(x) cos(5?/) dominates 
the large scale flow. The evolution of the energy of the large scale flow that is shown 



f Note that here refers to the critical energy input rate when hyperdiffusion is taken into 
account, which is four times greater than the critical input rate with u = 0. 
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Figure 6. The same as in figure 5, but for /3 = 10. 



in figure 9(b) increases rapidly and eventually saturates after about t = 150/cr(i 5). At 
this point the large scale flow gets attracted to a translating finite amplitude equilib- 
rium structure close in form to the harmonic Z = cos(x) cos(5y) (cf. figure 9(c)), drifting 
westward. The HovmoUer diagram of ip{x, y = tt/4, t), shown in 9(d), illustrates that the 
phase speed of the translating equilibrium structure is approximately equal to the phase 
speed of the unstable (n, m) = (1, 5) eigenmode that is also indicated in the figure. 

Consider now the energy input rate e = lOec for which the growth rates are shown in 
figure 8(b). While the maximum growth rate still occurs for the (|n|, |m|) = (1,5) non- 
zonal structure, zonal jet perturbations are unstable as well. In previous studies of SSST 
dynamics restricted to the interaction between zonal flows and turbulence, these initially 
unstable stationary structures were found to equilibrate at finite amplitude (Srinivasan & 
Young 2012; Constantinou et al. 2013). However, in the context of the generalized SSST 
analysis in this work that takes into account the dynamics of the interaction between 
coherent non-zonal structures and jets, we find that these SSST jet equilibria can be 
saddles: stable to zonal jet perturbations but unstable to non-zonal perturbations. To 
show this, we consider the evolution of a small zonal jet perturbation 6Z = cos(5?/) that is 
shown in figure 10. The initial zonal jet perturbation grows exponentially and its energy 
saturates at about t = 500 (a snapshot of the streamfunction at this time is shown at the 
left inset in figure 10). But soon after, non-zonal undulations appear and start to grow 
and the flow transitions to the stable Z = cos(a;) cos(5y) translating equilibrium state 
that is also shown in figure 10. As a result, the finite equilibrium zonal jet structure is 
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Figure 7. Dispersion relation of the unstable modes for /3 — 100. Growth rate o"r as a function 
of wavenumbers {n,rh) at (a) e = 2ec and (b) e = 50ec. Only positive values are shown. Ratio 
of the frequency of the unstable modes Gi over the corresponding frequency of a Rossby wave 
with the same wavenumbers an at (c) e = 2lc and (d) e = 50ec. Values of one denote an exact 
match with the Rossby wave phase speed. 
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Figure 8. The growth rate, Or as a function of the wavenumbers |n|, \m\ of the coherent 
structures at (a) e/tc = 4 and (b) e/tc = 10 (only positive values of Or are shown). 



structurally unstable to coherent non-zonal perturbations and is not expected to appear 
in the non-hnear simulations despite the fact that the zero flow equilibrium is unstable 
to zonal jet perturbations. We will elaborate more on this issue in the next section. 

Finally, consider the case e — SOec- At this energy input rate, the fast growing non-zonal 
perturbations cannot equilibrate, as the finite amplitude non-zonal translating equilib- 
rium states become SSST unstable. The unstable zonal jet structures on the other hand 
do equilibrate but their structure at finite amplitude is not zonally symmetric. The equi- 
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Figure 9. Equilibration of the SSST instabilities, (a) Streamfunction of the initial perturba- 
tion, (b) Energy evolution of the initial perturbation shown in panel (a) as obtained from the 
integration of the SSST equations (3.5) and (3.10) (dashed line) and from the integration of 
the quasi-linear system (3.1)-(3.2) with Nena ~ 10 (solid hue) and Nena ~ 100 (dash-dotted 
line) ensemble members, (c) Snapshot of the streamfunction v&eq of the translating equilibrium 
structure and (d) HovmoUer diagram of '^^q(x,y = 7r/4,f) for the equilibrated structure. The 
thick dashed line shows the phase speed obtained from the stability equation (4.2). The energy 
input rate is e = 4ec and /3 = 100. 



librium with the largest domain of attraction that is shown in figure 11, consists of a 
large amplitude zonally symmetric jet with small amplitude non-zonal propagating vor- 
tices embedded in it. These vortices that are shown in figure 11(b) to have approximately 
the compact support structure ^' = cos(2a;) cos(6j/) propagate westward as revealed by 
the HovmoUer diagram of the streamfunction shown in 11(c). 



6. Comparison to non-linear simulations 

The autonomous non-linear system of the SSST and its fixed points formally exist 
only in the infinite ensemble limit, in which the fiuctuations induced by the stochastic 
forcing are averaged to zero. In realistic atmospheres, such an approximation may not 
be valid. However, the characteristics of the stable equilibria as well as the occurrence 
of structural instability manifests even in single realizations of the turbulent system. For 
example, previous studies using SSST obtained zonal jet equilibria in barotropic, shallow 
water and baroclinic flows in close correspondence with observed jets in planetary flows 
(Farrcll & loannou 2007, 2008, 2009c, a). In addition, previous studies of SSST dynamics 
restricted to the interaction between zonal flows and turbulence in a /3-plane channel 
showed that when the energy input rate is such that the zero mean flow equilibrium is 
unstable, zonal jets appear in the non-linear simulations as well with the jet structure 
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Figure 10. Energy evolution of an initial jet perturbation SZ — 0.001 cos(53/) for e = lOtc and 
/3 = 100. The insets show a snapshot of the mean flow streamfunction ai t = 500 (left) and the 
streamfunction of the equilibrated structure at t = 6500 (right). 

(scale and amplitude) being accurately predicted by the theory (Srinivasan & Young 
2012; Constantinou et al. 2013; Tobias & Marston 2013). 

A very useful intermediate model that retains the wave-mean flow dynamics of the 
SSST system while relaxing the infinite ensemble approximation can be constructed by 
ignoring in (3.1) the non- linear term /"'. Then (3.1)-(3.2) become a quasi-linear system 
in which the ensemble mean can be taken over a finite number of ensemble members. Its 
integration is done as follows. Using the same pseudo-spectral code as in the non-linear 
simulations, A^ens separate integrations of (3.1) are performed at each time step with the 
eddies evolving according to the instantaneous flow. Then the ensemble average vorticity 
flux divergence is calculated and (3.2) is stepped forward in time according to those fluxes. 
The SSST equilibria manifest in the quasi-linear integrations with the addition of some 
'thermal noise' due to the stochasticity of the forcing that is retained in this system. This 
is shown in figure 9b where the energy growth of the coherent structure for Nens = 10 and 
Nens ~ 100 is plotted for the same parameters as the SSST integration. We observe that 
the energy of the coherent structure in the quasi-linear integrations fluctuates around 
the values predicted by the SSST system with the fluctuations decreasing as the square 
root of the number of ensembles. However, even with only 10 ensemble members we 
get an estimate that is very close to the theoretical estimate of the inflnite ensemble 
members obtained from the SSST integration. The translating equilibrium structure and 
its phase speed in the quasi-linear integrations are also in very good agreement with 
the corresponding structure obtained from the SSST integration (not shown). Since the 
quasi-linear system accurately captures the characteristics of the emerging structures 
yet is much faster to integrate than the corresponding SSST system, we will use the 
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(a) 




Figure 11. Equilibrium state with the largest domain of attraction for e = 40ec and /3 = 100. (a) 
Snapshot of the streamfunction of the equilibrium state, (b) Contour plot of the non-zonal 
component "fe, — 'i/^q of the equilibrium structure, where the overbar denotes a zonal average, 
(c) HovmoUer diagram of 9eq{x,y = n/4,t) for the equilibrated structure. 

quasi-linear integrations to test the predictions of the generalized SSST theory for the 
emergence and equilibration of both zonal and non-zonal coherent structures against the 
results of non-linear simulations. We will therefore integrate the quasi-linear system with 
Nens = 10 for the same parameters as in the non-linear simulations and compare the 
results. 

For the parameters chosen (/3 = 100), SSST predicts the emergence of non-zonal struc- 
tures when the energy input rate exceeds the critical threshold ec and the equilibration of 
the incipient instabilities into finite amplitude structures that should manifest in the non- 
linear simulations. The rapid increase of the nzmf index in the non-linear and quasi-linear 
simulations for e > ec shown in figure 1, illustrates that this regime transition in the flow 
is accurately predicted by SSST and that the quasi-linear and non-linear dynamics share 
the same bifurcation structure. While all stable SSST equilibria are in principle viable 
repositories of energy in the turbulent flow and the non-linear system is expected to visit 
their attractors for finite time intervals, significant power is expected to accumulate at the 
structure corresponding to the SSST equilibrium with the largest domain of attraction. 
Indeed for e = Ate the translating equilibrium with (|fc|, \l\) = (1,5) that has the largest 
domain of attraction, is the dominant structure in the non-linear simulations. Compar- 
ison of the energy spectra obtained from the quasi-linear and the non-linear simulation 
shown in figures 12a and 2a respectively, reveals that the amplitude of this structure in 
the quasi-linear and in the non-linear dynamics almost matches. Remarkably, the phase 
speed of the SSST translating equilibrium matches as well with the corresponding phase 
speed of the (|/c|, \l\) = (1,5) structure observed in the non-linear simulations, as can be 
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(a) (b) 




Figure 12. Time averaged energy power spectra, log(_E(fc, /)), obtained from the quasi-linear 
simulation of (2.1) at (a) e/ec = 4 and (b) e/ec = 50 or e = 3.3eni. There is a very good agree- 
ment for both energy input rates with the corresponding spectra obtained from the non-linear 
simulation (cf. figure 2). 

seen in the HovmoUer diagram in figure 3(b). Such an agreement in the characteristics 
of the emerging structures between the quasi-linear and non-linear simulations occurs 
for a wide range of energy input rates as can be seen by comparing the nzmf indices in 
figure 1. As a result, SSST can accurately predict the dominant non-zonal propagating 
structures in the non-linear simulations, as well as their amplitude and phase speed. 

The second transition in the flow that occurs with the emergence of zonal jets is 
more intriguing. The stability equation (4.2) predicts that the zonal structures become 
unstable at esz — 5.2ec. As discussed in the previous section however, the unstable zonal 
jets cannot equilibrate and the flow transitions to a non-zonal translating equilibrium 
state (cf figure 10). We found that structural instability of the finite amplitude zonal 
jet equilibria occurs for energy input rates in the range < e < e„;. When e > e„i, 
the non-zonal translating equilibrium states become SSST unstable (Constantinou et al. 
2013) while at the same time the zonal jet equilibria similar to the one shown in figure 
11 are stable and zonal jets emerge in both non- linear and quasi-linear simulations. The 
power spectra obtained from the quasi-linear simulations for e = 3.3e„; shows that both 
the scale and the amplitude of the zonal jets is accurately captured by SSST. Such 
an agreement again holds for a wide range of energy input rates, as the zmf indices 
obtained from the quasi-linear and the non-linear simulations indicate. Therefore, SSST 
can accurately predict the characteristics of both non-zonal propagating structures and 
of zonal jets in the non-linear simulations. 

7. Conclusions 

The emergence of zonal jets and non-zonal coherent structures out of a homogeneous 
background of turbulence in a barotropic beta-plane is investigated in this work. Non- 
linear simulations of a non-divergent fiow in a doubly periodic channel subjected to 
a homogeneous and isotropic random stirring that injects energy to a narrow band of 
wavenumbers, show two major flow transitions as the energy input rate increases. In the 
first, the translational symmetry in the flow is broken with the emergence of propagating 
coherent non-zonal waves that approximately follow the Rossby wave dispersion. The 
power in these non-zonal structures increases with the energy input rate until the second 
transition occurs with the emergence of robust zonal jets. While after the second tran- 
sition the zonal jets contain over half the energy in the fiow, significant power remains 
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in the non-zonal structures that remain coherent but propagate slower in the retrograde 
direction compared to Rossby waves. 

The two flow transitions and the characteristics of both the non-zonal structures and 
the zonal jets are then investigated using the tools of Stochastic Structural Stability The- 
ory (SSST). In the context of SSST, the cumulants of the dynamical variables containing 
the information on the ensemble mean values of the variables and their correlations are 
formed. The turbulent flow dynamics are then expressed as a systematic expansion in 
the cumulants that is truncated to second order. This truncation is equivalent to the 
quasi-linear dynamics formed by retaining the interaction between the eddies with the 
instantaneous mean flow, while neglecting the fast eddy-eddy interactions. With the in- 
terpretation of the ensemble average as a Reynolds average over the fast turbulent eddies 
adopted in this work, the second order cumulant expansion results in a closed, non- 
linear dynamical system the governs the joint evolution of the slowly varying coherent 
structures and the second order statistics of the rapidly evolving turbulent eddies. 

The stability of the fixed point of the non-linear SSST system with homogeneous 
enstrophy covariance and no mean velocity was examined. Structural instability was 
found to occur for perturbations with smaller scale than the forcing, when the energy 
input rate e = eKj/r^ is larger than a certain threshold ic that depends on /3 = /3/rKf. 

It was found that when /3 is small or order one, both zonal jets and non-zonal structures 
are unstable when the energy input rate is larger than ic, with the maximum growth 
rate occurring for stationary zonal structures. When /? is large, non-zonal structures first 
become unstable as the input rate increases past Sc with zonal jet structures becoming 
unstable at larger energy input rates. The maximum growth rate occurs in this case for 
non-zonal structures that propagate in the retrograde direction. These waves follow the 
Rossby wave dispersion for low supercritical values of the energy input rate, but propagate 
with significantly lower phase speeds at higher supcrcriticality. The equilibration of the 
unstable, exponentially growing coherent structures for large /3 was then studied. When 
the forcing amplitude is slightly supercritical, the finite amplitude equilibrium with the 
largest domain of attraction has a structure close to the corresponding unstable non-zonal 
perturbation with the same scale. 

The predictions of SSST where then compared to the results from the non-linear sim- 
ulations. The critical threshold above which coherent non-zonal structures are unstable 
according to the stability analysis of the SSST system was found to be in excellent agree- 
ment with the critical value above which non-zonal structures acquire significant power 
in the non-linear simulations. The scale, phase speed and amplitude of the dominant 
structures in the non-linear simulations were also found to be accurately predicted by 
SSST as the scale, phase speed and amplitude of the translating SSST equilibrium with 
the largest domain of attraction. In addition, the threshold for the emergence of jets 
that is identified in SSST as the energy input rate at which an SSST stable, finite am- 
plitude zonal jet equilibrium exists, was found to match the corresponding threshold for 
jet formation in the non-linear simulations, with the emerging jet scale and amplitude 
being accurately obtained using SSST. In summary, it was found that SSST predicts 
the two regime transitions in the turbulent flow as the energy input rate is increased: 
the symmetry breaking of the homogeneous equilibrium with the emergence of coherent, 
propagating non-zonal structures and the emergence of zonal jets. It also predicts the 
characteristics of the emerging structures (their scales and their phase speed), as well as 
their amplitude. 
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Appendix A. Calculation of the dispersion relation and its properties 

In this Appendix we derive the dispersion relation (4.2), which determines the stability 
of zonal as well as non-zonal perturbations in homogeneous turbulence. We follow closely 
the treatment of Srinivasan & Young (2012). We first rewrite (3.5), (3.10) in terms of 
the variables x = xi — X2, x = (l/2)(xi + X2), y = yi — 1/2 and y ~ {l/2){yi + 7/2). The 
derivatives transform into this new system of coordinates to = {l/2)dx + (— l)'^^i9i. 



% = (1/2)%- 



dh + and A = 9|j + 



(-1)^+155, A, = A + (1/4)A + i-iy+'dl^ + (-1) 



'+^d'i^, with A 



covariance S{x, x, y, y) 



It is also convenient to introduce the streamfunction 
{fp'ii}2), which is related to C{x,x,y,y) via: 



1- 



yyj 



S. (A 1) 



Equations (3.5), (3.10) then become in the absence of hyper-viscosity {ly = 0): 



dt + Uds + Udi + Vdy + Vdy 



C 



2(/3 + Zy)di -I- i^a^ - 2Z,% - 



(/3 + Zy)d^ + Zydi - Zxdy - Z.^dy 
(a|^ + ay 5=-2rC + S, 



A + -A 1 5 



(A 2) 



dtZ + UZ, + V{f5 + Zy) = idly - dl^){dl^ + dly)SU^y=o - rZ, (A 3) 

where {U ,V) - (l/2)(;7i + U2,Vi +¥2), {U ,V) = (C/i - C/2,Fi - V2), (Z^.Zy) = 
(1 /2){Zxi + Zx2 , Zy^ + Zy^ ) and {Zx, Zy) = (Zx-^ — Zx2: Zy^ — Zy^). 

The forcing covariance S is homogeneous and as a result it depends only on the dif- 
ference coordinates, x and y. It can then be readily shown from (A2)-(A3), that the 
state with no coherent flow [U^ = = Z^ = 0) and with the homogeneous vortic- 
ity covariance C^{x,y) = S/2r (implying also that the streamfunction covariance S is 
homogenous) is a fixed point of the SSST system. The stability of this homogeneous 
equilibrium, can be addressed by first linearizing the SSST system about it: 



dt6C 



dtSZ 



SUdx + 6Vdy ) C 
1 



SZyds-SZxdy) AS 



rA 



ckr-2{dl^ + dly)d,}SS^2r6C\ 



(A 4) 
(A 5) 



+ (dljj ~ 5|^)(9L + dySSU^y=o - rSZ, 

where SZ, 511, dV, SZx, SZy, SC and 6S are small perturbations in the ensemble mean 
vorticity, velocities and vorticity gradients and in the eddy vorticity and streamfunction 
covariances respectively, and then performing an eigenanalysis of the linearized equations 
(A4)-(A5). 

We consider a harmonic vorticity perturbation of the form SZ = e'"^+'™J'e°'* , for which: 



[6U,dV,dZx,6Zy] = -2 



(A 6) 
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Taking the same form for the streamfunctfon covariance perturbatfon 6S = Snm{x, y)e'"^+""^e'^* 
and inserting it in (A4)-(A5) along with (A 6) yields: 



(ct + 2r) 



A 

4 



nx 



my 



{mdi - ndy) 



- 2il3A+di + i7i/3 ( A - — 



AT^A 



(cr + r)N'^ + mfi = {mdi - ndy) A+S'„™|s=^=o, 



(A 7) 



(A 8) 



nidy and — - ' 



where N = n + m , A+ — nd- 
vorticity covariance with zero mean flow. 

Defining the Fourier transform of Snmix, y) by 



^/2r = A^S"^ is the equilibrium 



1 

2^ 



OO nOO 



-ikx — i 



^ydxdy 



(A 9) 



— CO ^ — OO 



we obtain from (A 7) that the Fourier component Snm satisfies: 

(mfc_ - nl_)Kl{Kl/N'^ - l)5'^(fc_,/_) 



>5'ti. 



2iPk{kn + ml) - m/3{Kl + Kl)/2 + {a + 2r)K\K': 



{mk+-nl+)Kl{KllN^ - l)S'^{k+,l+) 



2iPk{kn + ml) - inP{Kl + K'^_)/2 + (cr + 2r)K\K': 



(A 10) 

S/[2rA'4] 



with fc± = A: ± n/2, /± = / ± m/2, = 4 + Z| and K'^ ^ k^ + P. 
is the Fourier transform of S^, and S is the Fourier transform of S. In addition, (A 8) 
becomes: 



pr-Z fOO i-OO 

inP- ((T + r)N^ = / / [nmik"^ - l^) + {m^ -n'^)kl] Snmdkdl = A+- A., 

(All) 

where 

_ ^ n f°° [nmik^ - l^) + - {mk± - nl±)Kl{Kl - N'^)S'^ {k±,l±) ^^^^ 

^^2^ 2i/3fc(fcn + ml) - ml3{Kl + Kl)/2 + (<t + 2r)KlKl 

(A 12) 

Equation (All) can be further simplified by noting that because the choice of xi and 
X2 is arbitrary, the forcing covariance satisfies the exchange symmetry S(a;i, X2,yi, 2/2) = 
S(x2, 3:1,2/2, 2/1 )■ In terms of the new variables, the exchange symmetry is written as 
'E.{x ,x , y ,y) = 2(— — y, y), and consequently S satisfies S(— fc,— Z) = S{k,l). As a 
result: 

A+ = -A_. (A 13) 

Using (A 13) and shifting the axes in the resulting integrals (fc — > fc+n/2 and / — > l+m/2), 
reduces (All) to the following dispersion relation: 

°° r°° (mk - nl) [nm{kl - l\) + (m^ - n'^)k+l+\ K^{K^ - N^)S^{k, I) 
-00 J -00 2if3k+{k+n + l+m) - in/3 {K^ + K^) /2+{a + 2r)K^K^ '^'"^^ 
7r(cr + r)Ar2 - i7rn/3, (A 14) 

where — {k + n)^ + {l + m)'^. The corresponding dispersion relation on a periodic box. 
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can be readily calculated by simply substituting the integrals in (A 14) by finite sums of 
integer wavenumbers. 
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